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We use the axially-defomied Skyrme Quasiparticle Random Phase Approximation (QRPA) to- 
gether with the SkM* energy-density functional, both as originally presented and with the time-odd 
part adjusted to reproduce the Gamow- Teller resonance energy in '^''^Pb, to calculate the matrix 
elements governing the neutrinoless double-beta decay of ^^Ge, ^^°Te, ^^^Xe, and ^*°Nd. Our ma- 
trix elements in ^'^"Te and ^^®Xe are significantly smaller than those of previous QRPA calculations, 
primarily because of the difference in pairing or deformation between the initial and final nuclei. In 
■"^Ge and ^'^"Nd our results are similar to those of less computationally intensive QRPA calculations. 
We suspect the ^®Ge result, however, because we are forced to use a spherical ground-state, even 
though the HFB indicates a deformed minimum. 



PACS numbers: 21.60.Jz, 23.40.Hc 

I. INTRODUCTION 

Neutrinoless (Oj^/3/3) double-beta decay can occur if 
neutrinos are Majorana particles, at a rate that depends 
on a weighted average of neutrino masses (see Refs. [Tl[5] 
for reviews). The experimental search Oi>/3p is approach- 
ing sensitivity to neutrino masses below 100 eV. Extract- 
ing a mass from the results, however, or setting a reliable 
upper limit, will require accurate values of the nuclear 
matrix elements governing the decay, matrix elements 
that cannot be measured and must therefore be calcu- 
lated. A number of theorists have attempted the cal- 
culations, applying several distinct methods. Among the 
most popular is the proton-neutron quasiparticle random 
phase approximation (QRPA). 

The QRPA can be carried out at various levels of so- 
phistication. So far, with only a few exceptions O 0], 
the mean-fields on which the QRPA is based have been 
spherical by fiat. In addition, they have never been con- 
sistent with the residual QRPA interaction, nor have all 
the nucleons ever been active degrees of freedom, free 
from confinement in an artificially inert core. Finally, 
even the active nucleons have been forced to occupy a 
few harmonic-oscillator shells and thus have never tasted 
the continuum. Here we overcome all these limitations, 
allowing axially symmetric deformation, using a modern 
and well-tested Skyrme functional for both the Hartree- 
Fock-Bogoliubov (HFB) mean-field calculation and the 
QRPA that is based on it, keeping all the nucleons ac- 
tive, and placing the nucleus inside a large cylindrical 
box, so that discretized versions of continuum states up 
to high energy are available. 

Deformed Skyrme-QRPA calculations of this type have 



been applied extensively in recent years to nuclear vi- 
brations (see e.g., [SHS]) and will soon be applied to 
single-beta decay [TU]. Our implementation, described 
in detail below, is via a B-spline-based HFB code with 
the above-mentioned cylindrical-box boundary condi- 
tions followed by the construction and diagonalization 
of the QRPA Hamiltonian matrix in the basis of canon- 
ical two-quasiparticle states. The calculations consume 
enough CPU hours to require a supercomputer, and so 
we restrict ourselves here to four isotopes — ^^Ge, ^^°Te, 
^^^Xe and ^^''Nd — used in the some of the most promis- 
ing of current or proposed experiments. The deformation 
and pairing in the initial and final nuclei are often quite 
different and matrix elements can be suppressed as a re- 
sult [3|; our numbers depend crucially on the overlap of 
intermediate-nucleus states created by exciting the ini- 
tial ground state with those created by exciting the final 
ground state. The QRPA supplies only transition am- 
plitudes and so must be extended to obtain the overlap. 
Here we will apply a prescription like that in Ref. [3], 
while noting that a well justified and tractable expres- 
sion is still lacking. 

This article is organized as follows: Section |TT] con- 
tains a brief overview of the matrix elements governing 
double-beta decay and of the Skyrme QRPA. Section [TTT| 
describes the details of our computational implementa- 
tion and Sec. |IV| presents our results. Section [V| is a 
conclusion. 



II. DOUBLE-BETA DECAY AND THE QRPA 
A. Decay operators 



The lifetime for 01/13/3 decay, if there are no heavy par- 
ticles mediating the decay, is 
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where (m^) is a weighted average of three neutrino 
masses, G'^'^ is a phase space factor (recently reconiputed 
in Ref. [H]), and M'^'^ is a nuclear matrix elemeniR Al- 
though the matrix element contains intermediate states 
and an energy denominator, it can to good approxima- 
tion |12| be represented by one involving only the initial 
and final ground states. In this "closure" approximation 
and neglecting the small tensor term, one can write the 
matrix element as 
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where the factor 1.25 is inserted by convention. Tab = 
l^a — ^b| is the distance between nucleons a and b, jo is the 
usual spherical Bessel function, E is an average excitation 
energy to which the matrix element is insensitive (and for 
which we use the value 10 MeV), and the nuclear radius 
R = 1.2A^/^ fm is inserted with a compensating factor 
in the phase-space function to make the matrix element 
dimensionless. The "form factors" hp and her are given 
by 
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Here rup and ttt-tt are the proton and pion masses. 

The two-neutrino double-beta decay {2iy/3/3) rate, 
which we will use to fit parameters for our cal- 
culation, can be written as 
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where G^" is another phase-space factor (also recom- 
puted in Ref. ^11^) and M'^" is a matrix element. The 



^ This matrix element differs from the unprimed M*^" used else- 
where by a factor of The two are equivalent when 
gA is taken to be 1.25, but differ when it is modified. (Actu- 
ally, qa is closer to 1.27 than 1.25, but we follow tradition here.) 
The convention we use puts all the qa dependence in the matrix 
element and none in the phase-space factor. 



closure approximation is not good for two-neutrino de- 
cay, and the matrix element must contain intermediate 
states explicitly: 
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where n labels states in the intermediate nucleus with 
energy E^ Mi and Mf are the masses of initial and fi- 
nal nuclei, and the effects we've neglected — forbidden 
currents, the Fermi matrix element, etc. — are small. 

Recent study [13l [M] has shown that realistic short- 
range correlations have only a small effect on the double- 
beta matrix elements. Including them here, even approx- 
imately, would complicate our computational procedure 
considerably and so we omit them altogether. 



B. Deformed Charge-changing QRPA 

The self-consistent axially-symmetric Skyrme-HFB- 
QRPA method for like-particle excitations, on which our 
code is based, is described thoroughly in Refs. [TS], [S], 
and 0. We modify the code discussed there in a rather 
straightforward way — changing the basis of like-two- 
quasiparticle states to a basis of one-quasiproton-one- 
quasineutron states, removing the Coulomb interaction, 
and keeping only the relevant parts of the Skyrmc func- 
tional — to work with charge-changing modes rather 
than like-particle modes. 

We adopt the Skyrme functional (or effective interac- 
tion) SkM* [16]; that functional has been shown to de- 
scribe nuclear deformation well and reproduces low-lying 
quadrupole vibrations in rare-earth nuclei noticeably bet- 
ter than the comparably popular functional SLy4 [S]. 
We modify the time-odd particle-hole part of the func- 
tional as in |17j . which discussed charge-changing tran- 
sitions, by setting the parameters (defined in that refer- 
ence) = 0, CY' = 0, and Cf [0] = Cf [p„m] = 100 Mev 
fm'^ (pnm is nuclear-matter density). With these modi- 
fications, the functional reproduces [TU] the location of 
the Gamow- Teller resonance and the fraction of observ- 
able strength in the resonance. We will report results 
with and without the modifications to show their effect. 

For the particle-particle part of the functional we use 
a simple volume (zero-range) pairing interaction, the 
strength of which we adjust separately in the isoscalar 
channel (T = 0) and in each of the three isovector (T = 1 
with Tz ~ —1, 0, and 1) channels. We describe the ad- 
justment in more detail in the next section. 

Evaluating the Oz//3/3 matrix elements requires a mul- 
tipole decomposition of M"^" suitable for cylindrical ge- 
ometry. The details of that appear in the Appendix. 



III. COMPUTATIONAL IMPLEMENTATION 

Only recently have fully self-consistent deformed 
Skyrme-QRPA calculations entered the scene. The com- 
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bination of methods we use here requires many thousands 
of CPU hours. Our methodology wiU at some point be 
obsolete because of the development of much faster Fi- 
nite Amplitude JTE' and iterative Arnoldi approaches, 
which use mean-field codes with time-independent con- 
straints to solve the QRPA equations. Our method, by 
contrast, involves the explicit construction and diagonal- 
ization of the QRPA Hamiltonian matrix in a basis of 
two-canonical-quasiparticle states. These states are ob- 
tained from the HFB calculation mentioned previously. 

To solve the initial HFB equations we use the Van- 
derbilt deformed HFB code [20], which represents wave 
functions in a basis of B-splines. Our cylindrical box has 
dimensions fmax = -^max — 20 fm, about three times as 
large as the radius of the heaviest nucleus studied here 
and a number found suitable in Ref. [20] . Our mesh 
spacing is 0.7 fm and the energy cutoff of the HFB solu- 
tions is 60 MeV. We do not restrict the deformation (ex- 
cept to be axially symmetric) but rather allow the mean 
field to evolve freely to the nearest local binding-energy 
minimum. Using a range of quadrupole deformation pa- 
rameters /32 as initial guesses, we find one or more local 
minima and select the most bound solution as the mean 
field on which we base the QRPA. In ''^Ge, however, we 
do not use the most bound solution; we discuss the rea- 
sons for this exception in the next section. To obtain 
the strength of the proton-proton and neutron-neutron 
{T = l.Tz = ±1) pairing interaction we match the HFB 
pairing gaps with the experimental pairing gaps obtained 
from a three-point interpolation formula, with separation 
energies from the ENSDF [21] database. 

The computational requirements for running our 
charge- changing QRPA code are significantly less than 
those for the like-particle code on which it is based be- 
cause a) the proton-neutron two-quasiparticle basis is 
only about half the size of the like-two-quasiparticle ba- 
sis, and b) the removal of the Coulomb interaction re- 
lieves us of a large computational burden. For a given 
multipole, our charge-changing code typically runs much 
faster than the like-particle code. That speedup, how- 
ever, still leaves us with runs that consume many thou- 
sands of CPU hours per multipole in each nucleus. 

We cannot include all one-quasiproton-one- 
quasineutron states in our QRPA basis and so truncate 
the same way as in Ref. [5] . The truncation is controlled 
by two parameters v^^^ and v^^^^. The first allows us to 
remove two-quasiparticle states with almost completely 
two-particle or two-hole nature (i.e. states that primarily 
lie in the AzL2 neighbors of the reference nucleus instead 
of the in intermediate double-beta nucleus), and the 
second lets us cut out states in which one of the particles 
is far below the Fermi surface and the other far above 
it. Such excitations have very high energy and do not 
mix significantly with lower-energy states. In practice 
15,000 two-quasiparticle states for the lowest multipoles, 
out of a total of about half a million, are enough to 
approximate the exact answer very well, making the 
construction and diagonalization of the QRPA matrix 



tractable on a supercomputer. 

After diagonalizing the QRPA Hamiltonian, we need 
to determine the double-beta-decay matrix elements. For 
QvjSp decay, the matrix element can be written as 
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where c|, are particle-creation operators, the indices with 
p refer to protons and those with n to neutrons, each 
index stands for the set of quantum numbers p — 
{jp,T^p,kp} (angular- momentum along on the intrinsic 
axis, parity, and an additional enumerating index), a mi- 
nus sign in front of an index means that the sign of the 
quantum number is reversed, and 
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The two-particle states in Eq. ([8]) are antisymmetrized. 
We use a multipole expansion, detailed in the Appendix, 
to evaluate the two-body matrix elements in Eq. ([s]) with 
B-spline integration. The coding for the two-neutrino 
two-body matrix elements, which we use to evaluate the 
matrix element in Eq. Q, requires no Bessel function 
expansion. 

Two-neutrino decay is simpler for another reason as 
well; only states with angular momentum and parity 
J'^ = 1"*" contribute to the matrix element. In our de- 
formed calculation we follow the usual procedure of rep- 
resenting laboratory states in a rigid-rotor approximation 
as combinations of a) Wigner functions Dfjj^ and D'lj_j^ 
of Euler angles, and b) an intrinsic QRPA state with a 
well-defined projection K along the symmetry axis of the 
angular momentum J. M'^'^ thus gets contributions only 
from states with \K\ < 1. In neutrinoless decay, on the 
other hand, states with any K'^ contribute. The contri- 
butions get progressively smaller as K gets larger. In- 
cluding state with \K\ < 10 is enough to approximate 
the matrix element accurately, as Fig. [l] shows. 

One interesting feature of Eq. Q is the presence of 
the overlap {N\N'). The QRPA is a small-amplitude ap- 
proximation and although it provides transition densities 
from a ground state to excited states, it can't, without 
extension, provide excited state wave functions. The ex- 
cited states \N) and |A^') are based on different quasi- 
particle vacua and the quasiboson approximation that is 
inherent in the QRPA erases the information necessary 
to relate the two vacua. Two expressions for the overlap 
have been given in the past few years: one, from Ref. 
[3], neglects "scattering terms" even though they can- 
not be shown to be small and the other, laid out in Ref. 
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FIG. 1. (Color online.) The cumulative '^'^Ge Oz//?/? ma- 
trix element (using SkM* and gA = 1-0) as the number of 
intermediate-state multipoles K'" is increased. Convergence 
is reached by \K\ — 10. Both positive and negative parities 
are included, as are both Fermi and Gamow- Teller contribu- 
tions. 







Ref. [24] Exp. 




this work 


Sk3 


SG2 Ref. [25] Ref. [26] 


^«Ge 


0.18^ 


0.161 


0.157 0.095(30) 0.2623(9) 


^^Se 


-0.018 


-0.181 


-0.191 0.163(33) 0.3090(37) 


i30Te 


0.01 


-0.076 


-0.039 0.035(23) 0.1184(14) 


i3"Xe 


0.13 


0.108 


0.161 - 0.1837(49) 


i36Xe 


0.004 


0.001 


0.016 - 0.122(10) 


i36Ba 


-0.021 


0.009 


0.070 - 0.1258(12) 




0.27 


0.266 


0.271 0.367(86) 0.2853(21) 




0.22 


0.207 


0.203 0.230(30) 0.1931(21) 



-0.025 used 



TABLE I. The quadrupole deformations (32 of the initial and 
final nuclei in our work, compared with the values obtained 
in [2j and experimental values from |25l 126] 



IV. RESULTS AND DISCUSSION 



[55] , uses the form of the boson vacuum but replaces the 
bosons with the fermion pairs from which they stem. Un- 
fortunately, this last idea leads to expressions that can 
only be evaluated perturbatively; these become unwieldy 
after the lowest couple of orders in the expansion, the 
convergence of which may not be fast. Here we sim- 
ply evaluate the overlap in the quasi- Tamm-Dancoff ap- 
proximation (neglecting the QRPA "Y" amplitudes); in 
this limit of the QRPA, excited states are well-defined 
two-quasiparticle excitations of HFB vacua and require 
no bosonization. The results are not very different from 
those obtained in the scheme proposed in Ref. [3|. We 
provide more details in the Appendix. 

As is typical in QRPA calculations, we use the mea- 
sured values of two-neutrino decay rates to fit proton- 
neutron pairing strengths. Following a suggestion in Ref. 
[23] . we adjust the isovector (T = l,Tz — 0) strength 
so that the Fermi 2i>/3/3 matrix element vanishes, as it 
should (almost) because the ground state of the final 
nucleus has a different isospin than the double-isobar- 
analog state of the initial nucleus. If instead we fix 
the proton-neutron isovector pairing strength at the av- 
erage of the proton-proton and neutron-neutron pairing 
strengths, we find a nearly identical result. The isoscalar 
pairing strength, which we call Vp here, is the parameter 
typically called gpp in other QRPA calculations. We ad- 
just it so as to reproduce the experimental two-neutrino 
matrix element, with both an unquenched {gA = 1.25, 
see the footnote) and quenched {gA — 1.0) axial- vector 
coupling constant. We then use the resulting pairing 
strengths in computing the Ov/3/3 matrix elements, once 
for each value of gA- For ^^°Te and ^'^^Xe, we compute 
the neutrinoless double-beta-decay matrix element with 
the unmodified SkM* over a range of isoscalar pairing 
values Vq to assess its sensitivity to the fit. 



We start by comparing the quadrupole deformation pa- 
rameters f32 obtained from our HFB calculation to other 
theoretical and experimental values in Tab. |l] With the 
exception of ^^Se, where our Skyrme-HFB computation 
fails to converge to a prolate solution, our quadrupole 
deformations are similar to those obtained using Sk3 and 
SG2 Skyrme interactions in Ref. 24J. The failure to con- 
verge is most likely due to a very flat bottom of the bind- 
ing energy curve with respect to deformation in "^^Se. 

In ^^Ge, the minimum energy occurs at a prolate de- 
formation of /32 — 0.18. This deformation is so dif- 
ferent from that of ''^Se, however, that our predicted 
two-neutrino matrix element is smaller than the mea- 
sured value no matter what we use for gA or Vq. We 
therefore choose to use the local near-spherical mini- 
mum {(32 — —0.02 5) for ^^Ge instead. As we shall 
see, this gives us a result that is not too different from 
other QRPA numbers, including those of Ref. [4, which 
presents both spherical-spherical and prolate-prolate cal- 
culations. It also indicates, however, that the QRPA is 
inadequate in this system. The soft surfaces with multi- 
ple minima require a formulation that mixes mean fields, 
e.g. the generator-coordinate method (often referred to 
as energy-density functional (EDF) theory) of Ref. [27] . 
or an extension thereof. 

In the daughter nucleus ^^°Xe we get a prolate solution, 
making ours the first QRPA calculation to take the de- 
formation into account in the decay of ^^°Te. The study 
in Ref. ,28J , using HFB with the Gogny interaction, finds 
a second minimum with oblate deformation and a barrier 
of only 1 MeV or so separating the two minima. As in 
^^Se, therefore, the use of a single mean-field in the con- 
struction of the ^^°Xe ground state is somewhat suspect. 

We turn now to the matrix elements themselves. Fig- 
ure [2] displays the dependence of the 2^(3(3 matrix element 
on the isoscalar pairing strength Vq in the four systems we 
study. We use the recent evaluation of the phase-space 
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FIG. 2. (Color online.) The dependence of two-neutrino 
double-beta decay matrix elements on Vo, the isoscalar pair- 
ing strength. The thick solid and dashed (red) curves are pro- 
duced by the original SkM* interaction and the dotted and 
thin (blue) curves by the modified interaction. The thick solid 
and dotted curves are computed with qa = 1.25, the dashed 
and thin solid curves with the quenched value qa ~ 1.0. 



factors in Ref. to extract the experimental matrix 
elements. Because M'^'^ for ^'^^Xe was just measured for 
the first time by the EXO-200 [29^ and KamLAND-Zen 
[5Uj experiments, ours is the first QRPA double beta de- 
cay computation to use an experimentally obtained value 
rather than an upper limit to determine the strength of 
isoscalar pairing. 

Figure [3] illustrates the dependence of the Oi//3f3 decay 
matrix element on Vq . The neutrinoless matrix element is 
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FIG. 3. (Color online.) The dependence of M™"' in ""Xe and 
^^°Te on the Vo produced by the unmodified SkM* interaction. 
The solid curve represents the results with gA = 1.0 and the 
dashed curve represents the results with gA = 1.25. 





gA = 1.0 


gA = 1.25 


gA = 1.0 


<7A = 1.25 


^^Ge 


4.40 


5.53 


4.12 


5.09 




1.13 


1.38 


1.32 


1.37 


i3«Xe 


1.26 


1.68 


1.18 


1.55 


"•^Xe (HFLN) 


1.54 


2.05 


1.44 


1.89 




2.52 


3.14 


2.14 


2.71 



TABLE II. The Oj//3/3 matrix elements in our Skyrme-HFB- 
QRPA calculation, with both the functional SkM* and a mod- 
ified version of it, and with both a quenched and unquenched 
axial- vector coupling constant gA- The second row of ^'^^Xe 
numbers contains the results with the Hartree-Fock -f Lipkin- 
Nogami overlap (see text). 



less sensitive to this pairing mode than the two-neutrino 
matrix element. We collect our final results for the Qvfip 
matrix elements with both gA = 1.25 and gA = 1.0 in 
Table |llj The modification of SkM* usually suppresses 
the Oi//3/3 matrix element, by up to 15%. It actually 
seems to increase the matrix element in -'^■^^Te by 17% for 
gA = 1.0, but as Fig.[2]shows, the fitting procedure for Vq 
with gA = 1.0 gives an anomalously small value, and so 
that result must be taken with a grain of salt. In Table [Till 
we compare our values with the modified SkM* and qa = 
1.25 with earlier theoretical results. Our matrix elements 
for ''^Ge and ^^°Nd are in a good agreement with the 
spherical result for Ge and the deformed one for Nd in 
Ref. H]. For ^'^^Xe and ^^°Te we get noticeably smaller 
matrix elements than obtained in prior work, all of which 
was carried out in the spherical QRPA. Figure |4] displays 
the same information as the table graphically. 

The suppression we see in ^^°Te can be attributed to 
the deformation of the daughter nucleus. Previous QRPA 
calculations for ^^°Te [T31 [31] have assumed spherical 
symmetry. We've already mentioned, however, that a 
single minimum may not be adequate to represent the 
ground state of ^'^"Xe. We suspect that the complete 
neglect of deformation in previous work leads to a ma- 





present 


QRPA/T 


QRPA/J ISM IBM-2 PHFB EDF 


^"^Ge 


5.09 


5.30, 4.69* 


5.355 2.96 5.465 — 4.60 




1.37 


4.92 


4.221 2.81 4.059 4.66 5.13 


"«Xe 


1.89 


3.11 


2.802 2.32 2.220 — 4.20 




2.71 


3.34* 


— — 2.321 3.24 1.71 



TABLE III. Comparison of our Ovfifi matrix elements, 
from the modified SkM* functional and gA ~ 1.25, with 
those obtained from the interacting shell model (ISM, ^31 ), 
QRPA calculations by the Tiibingen group gUIS] (QRPA/t) 
and Jyvaskyla group [32] (QRPA/J), the energy-density- 
functional method (EDF), projected HFB [33] (PHFB) 
and the interacting boson model [34] (IBM-2). Prior results 
that include deformation are indicated by a star. In ^"^^Xe we 
use the Hartree-Fock + Lipkin-Nogami overlap to scale our 
matrix element. 
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FIG. 4. (Color online.) The results of Tablejlll : 
The solid (blue) arrow points to our result with the Hartree- 
Fock + Lipkin-Nogami overlap in ^'^^Xe. The dashed (red) 
arrow points to the new shell-model result of Ref. [35] in the 
same nucleus. 



trix element that is too large, but it may also be that 
our sharp prolate Xe ground state yields one that is too 
small. 

The other decay in which we disagree significantly with 
previous QRPA calculations is that of ^'^^Xe. Our signif- 
icantly smaller result here is not caused by deformation 
difference, nor does it come from the availability of new 
two-neutrino decay data. Instead, it can be traced to the 
overlap between the initial and final HFB mean fields. 
This overlap usually reflects the difference in deformation 
between the mother and daughter nuclei, and for that 
reason has been completely neglected in previous QRPA 
calculations for the decay of ^'^^Xe, where both the ini- 
tial and final nuclei are spherical. We find here, however, 
that differences in pairing structure in the neutron mean 
fields lead to a smaU overlap: (HFB/|HFB,) = 0.47. The 
suppression is related to the = 82 shell closure, which 
produces a sharp Fermi surface that smooths measurably 
with the addition of two neutrons. We see no reason to 
completely neglect the overlap, but the situation may be 
analogous to the that in the decay of ^^°Te. A more 
realistic representation of pairing than is offered by the 
HFB mean field might make the difference in structure 
between the initial and final nuclei a little less dramatic. 

To test that last assumption, we reevaluated the over- 
lap in the Hartree-Fock -I- Lipkin-Nogami approxima- 
tion [3^, which reduces fluctuations in particle number 
and prevents the total breakdown of pairing at closed 
shells. The new occupation numbers increase the overlap 



to 0.57. If we scale the matrix element with qa = 1-25 
and the modified Skyrme functional accordingly, we ob- 
tain the result 1.89, indicated by the arrow in Fig. 
|4] Although the Lipkin-Nogami calculation is not self- 
consistent — we have only modified the overlap, not the 
HFB quasiparticle and the QRPA strength calculations 
— this larger number is our best estimate of the ma- 
trix element. The deviation of the Lipkin-Nogami over- 
lap from unity, while less than that of the HFB overlap, 
still makes the result smaller than those of any other 
QRPA calculations. Interestingly, a recent shell-model 
[3 5) calculation finds that increasing the model-space size 
produces the smallest matrix element yet for this decay: 
1.46. 

All the substantial differences between our QRPA cal- 
culations and others can be traced to deformation or pair- 
ing effects that were neglected in previous work. Our use 
of a self-consistent QRPA with all nucleons treated as 
active participants, the continuum accounted for, etc., 
doesn't, in itself, change results dramatically. That find- 
ing is not altogether surprising. Self-consistency is im- 
portant in the QRPA partly because it eliminates spu- 
rious strength. In the charge changing QRPA, however, 
the absence of proton-neutron mixing in the HFB and 
the explicit breaking of isospin mean that there is no 
spurious strength even in non-self-consistent calculations. 
More importantly, it is already well known ^37j that dif- 
ferences between variants of the QRPA largely disappear 
when the strength of the isoscalar pairing interaction is 
adjusted so that each variant reproduces the measured 
two-neutrino rate. Our variant does not escape this fate; 
that one parameter is a like a broad and coarse brush that 
paints over any sophistication in the underlying method. 



V. CONCLUSIONS 

We have performed large-scale Skyrme-HFB-QRPA 
computations for four important double-beta emitters. 
We have allowed for axial deformation of the initial and 
final nuclei. Our implementation increases the scale of 
the computation to the limits of contemporary technol- 
ogy- 

For ^"^Ge and the very deformed ^^"Nd, our results are 
in line with the earlier results of Ref. [1]. We note, how- 
ever, that the assumption both here and elsewhere that 
the ^^Ge and ^^Se ground states are spherical probably 
results in a matrix element that is too large. 

In ^'^'^Te we improved on the ground state used in pre- 
vious QRPA calculations by taking into account the de- 
formation of the final nucleus. Shape coexistence in the 
daughter ^■^'^Xe is beyond the scope of the QRPA, how- 
ever; if present, it could further modify the value of the 
matrix element. 

Our ^'^^Xe matrix element is the first QRPA result ob- 
tained from the new two-neutrino-decay measurements. 
It is also the first to take into account the overlap of 
the two sets of QRPA intermediate states. The overlap 
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is smaller than one might expect because of the sharp 
neutron Fermi surface in the initial nucleus. In reality, 
the Fermi surface cannot be perfectly sharp, and the true 
matrix element is probably better represented by replac- 
ing the HFB overlap with the Hartree-Fock + Lipkin- 
Nogami overlap. The modified result is larger than our 
original matrix element but still smaller than those of 
other QRPA calculations. 

Our computation demonstrates that there is little to be 
gained by further increasing the size and sophistication 
of QRPA calculations. Any straightforward alterations 
to the QRPA, other than the development of a better 
energy-density functional, are unlikely to improve the re- 
sults substantially. We have reached the point at which 
shortcomings of the QRPA itself restrict improvement. 
The inability to treat shape coexistence is an issue at 
least for the daughter nuclei ^^Se and ^■^"Xe. The mean- 
field treatment of pairing may be a problem in nuclei such 
as ^^^Xe that have closed shells. We can address these 
issues only by moving beyond the QRPA. 
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Appendix: Neutrinoless double beta decay matrix 
elements in the cylindrical box 

Equation ([t]) is in essence a trace of a product of four 
large square matrices. The transition densities in that 
equation are 

1+1 J 



(Oj|cLpC„|7V) = SpVpUnXp„ + UpSnVnY_p_„ , (A.l) 



and 



(iV'|c|^,C_„'|0+) = -UpiSn'Vn'Xp,^, - Sp'Vp'Un'Yp 



N' 



rN' 



(A.2) 



Here the indices p and p' indicate protons, and n and 
n' neutrons, as discussed in the main text; c| is a pro- 
ton creation operator, and Up and SpVp are the proton 
occupation amplitudes in the canonical basis, in the no- 
tation of Ref . [38j . and Y^^_^ are forward-going and 
backward-going QRPA amplitudes. 

To reduce the number of nested numerical integrals in 
the Qv(3f5 matrix elements in Eq. ([T]), we take advantage 
of the following expansion for the spherical Bessel func- 
tion in Eq. (Isl): 



(qrab) ^ 4:Tr^ji{qra)jiiqrb) ^ Yi*^{fa)Yi,n{h) 



1=0 



rn= — l 



(A.3) 

This allows us to separate the integrals over coordinates 
of the two nucleons: 



dq 



(2/ + 1) 



x/i^„(g)/^f_„,(g), (A.4) 



and 



pn,p n 



dq 



qhQT{q^) 



1 



q 



E (2^ + 1) 



/— max(0,/C— /J.) 



(Z + (if-M))! 



X /; 



l,K-p.,-p,f \ jlM-p,p 



{q)lT-n'i^) 



(A.5) 



where Save — E — {Ei + Ef)/2. Naturally, the infinite 
summations over I must be truncated. For most values 
of the neutrino energy q, not many terms are needed 
for convergence. In the program we truncate the expan- 
sion dynamically by requiring a preset accuracy in the 
quadrature for each value of q. 

The axial symmetry of the normalized canonical single- 
particle wave functions means that they can be written 
in the form 



«'a(r) 



^ E V'a(s;p,z)e^(^'"-^)*x., (A.6) 



where s is the spin projection, Xs is a standard two- 
component spinor, and the angular-momentum pro- 
jection onto the intrinsic axis. The integrations over the 
azimuthal angle cj) are trivial and the integrals ij^iq) and 
l]^'^{q) are therefore only two-dimensional: 



oo /"OO 



labiq)^ dz dppi;l{p,z)MP,z) 



xji{q^^T^)PP 



(A.7) 
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and 



oo roo 



J -oo Jo 



(A., 



Here the P™(x) are the usual associated Legendre poly- 
nomials, a I, are the Pauli matrices in the spherical vector 
basis and 



(A.9) 



As discussed in the main body of the text, we also 
need to evaluate the overlaps {N\N') between the QRPA 
states stemming from different mean fields. A satisfac- 
tory expression for these is lacking, but the chief ingre- 
dient in any such expression will be the overlap of two 
HFB vacua. The generalized Thouless theorem [38^ re- 
lating the two non-orthogonal quasiparticle vacua |HFBi) 
(initial state) and |HFB/) (final state) to each other is: 



|HFB,) =AA-icxp 

\ kl 



„(/)t„(/)t 

kio-k ai 



I HFB/ 



(A.IO) 

where the a): are quasiparticle creation operators in the 
final nucleus. The normalization factor is related to the 



transformation coefficients D^i via the Onishi formula: 
TV = (HFB/|HFB,)-i y^dct(l + DW) . (A.ll) 

Because the canonical-basis wave functions form a com- 
plete set, there exists a linear transformation between the 
two HFB solutions: 



,(/)t _ 



where 



and 



= 5](7^,„a«t+5,._„a«) 



T^fen = {n\k){ukUn + SkVkSnVn) , 



Sk.-r 



{n\k) {uk 



(A.12) 



(A.13) 



(A.14) 



Substituting Eq. ( |A.fO ) and (A.12) into the definition of 
the quasiparticle vacuum 



aL^i|HFB;)=0, 



(A.15) 



expanding the exponential, and comparing the terms con- 
taining one quasiparticle creation operator, we get the 
matrix equation 



n*D 



-S* 



(A.16) 



from which we can obtain the transformation coefficients 
Dki. 

As mentioned earlier, we approximate the QRPA over- 
laps states by QTDA overlaps, i.e. by neglecting the F's. 
This leads finally to the expression 



{N\N') = AA-l ^ ^ X^:X^,'^, np,p + Sp'p"Dp„p n„'n + J2 ^n'n"D„„,, 



(A.17) 



pn p'n' 



This formula differs slightly from the one presented in 
Ref. [3] and used in most QRPA double-beta-decay calcu- 
lations. Our overlap differs in that we keep the transfor- 
mation between the two HFB bases accurate and neglect 
the usually tiny term proportional to two Y amplitudes. 
In test calculations we find the numerical difference be- 



tween the two prescriptions to be negligible, as the com- 
mon leading term is already a good approximation. A 
more consistent evaluation of these overlaps that includes 
ground-state correlations can easily get both very com- 
plicated and computationally demanding, as evidenced 
by recent work in the like-particle QRPA in Ref. |22) . 
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